---
title: "Replication File for Analyses in 'Political Knowledge in the Context of Changing Institutions'"
author: "Sofia Vidotto, Rebecca Weitz-Shapiro and Matthew S. Winters"
date: "`r Sys.Date()`"
output: html_document
---

# Note that Figures 1 and 2 in the manuscript, which summarize recent literature on political knowledge, are produced separately from this replication archive; here, we provide replication code for the substantive empirical analyses in the paper.

```{r}
library(tidyverse)
library(foreach)
library(survey)
library(haven)
library(dplyr)
library(plyr)
library(ggrepel)
library(readr)
library(ggplot2)
library(ggpubr)
library(cowplot)
library(gridGraphics)
library(grid)
library(gridExtra)
library(patchwork)
library(MASS)
library(sandwich)
library(lmtest)
library(lme4)
library(stargazer)
library(estimatr)
```

```{r}
#States LAPOP
plot_states_dt <- read_csv("Datasets/Processed datasets/lapop/states_means.csv") 
# Presidential term LAPOP
plot_term_lrc_dt <- read_csv("Datasets/Processed datasets/lapop/term_means.csv")
# No of Representatives LAPOP
plot_reps_right <- read_csv("Datasets/Processed datasets/lapop/reps_means.csv")
# Reelection/term limits Afrobarometer
reelection_dt <- read_csv("Datasets/Processed datasets/afrob/term_means.csv")

# Individual-level data 
dt2 <- read_csv("Datasets/Processed datasets/lapop/dt2.csv") 
dt <- read_csv("Datasets/Processed datasets/afrob/dt.csv")
dt3 <- read_csv("Datasets/Processed datasets/cses/dt3.csv")

# Read the files
party_dt <- read_csv("Datasets/Processed datasets/afrob/party_means.csv")
leader_dt <- read_csv("Datasets/Processed datasets/afrob/leader_means.csv")
cses_party_mean <- read_csv("Datasets/Processed datasets/cses/cses_party_mean.csv")
```

```{r}
######################
#### Figure 3 ########
######################

######################
## Number of states ##
######################

# Note on number of data points by country - for the question about the number of states, countries can have two, three or four data points. Countries with two data points correspond to surveys collected in 2008 and 2010 (Argentina, Bolivia, Brazil, and Paraguay), countries with three data points correspond to surveys collected in 2006, 2008 and 2010 (Chile, Dominican Republic, Haiti, Peru, Uruguay and Venezuela), and countries with four data points correspond to surveys collected in 2004,2006,2008, and 2010 (Colombia, Costa Rica, Ecuador, El Salvador, Guatemala, Honduras, Mexico, Nicaragua, and Panama)
  
#Calculate the mean
plot_states_dt <- 
  plot_states_dt %>%
  filter(!is.na(states_years_since_asinh)) %>%
  dplyr::select(pais, states, states_years_since, states_years_since_asinh, gdp_log, polity2, partipdem, delibdem, polyarchy)%>% 
  group_by(pais)%>%
  dplyr::summarize(states = mean(states, na.rm = T), states_years_since = mean(states_years_since,  na.rm = T), states_years_since_asinh = mean(states_years_since_asinh,na.rm = T), gdp_log = mean(gdp_log,na.rm = T), polity2 = mean(polity2,na.rm = T), partipdem= mean(partipdem, na.rm = T), delibdem= mean(delibdem, na.rm = T), polyarchy= mean(polyarchy, na.rm = T)) %>%
  ungroup() 

# Define the custom breaks 
custom_breaks <- c(0, 15, 45, 75, 105, 145) 

plot_states_dt %>%
mutate(
    abv = paste(toupper(substr(str_remove(pais, " "), start = 1, stop = 3)))
  ) %>%
  ggplot() +
  aes(
    x = states_years_since_asinh,
    y = states
  ) +
  geom_point(shape = 16) +
  geom_smooth(method = "lm") +
  geom_text_repel(aes(label=abv), size=2.5)+
  scale_x_continuous(breaks = asinh(custom_breaks), labels = custom_breaks) +
  scale_y_continuous(limits = c(0.2, 1), breaks = seq(0.2,1,by=0.2)) +
  labs(
    x = "Years since number of states changed",
    y = "Proportion of correct answers",
    title = "How many states does the country have?",
    caption = "Source: AmericasBarometer (2004-2014)"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot1_LA

print(plot1_LA)

#######################
## Presidential term ##
#######################

# Length of the presidential term
# Note on number of data points by country - For the question about the length of the presidential term, countries can have four, five, or six data points. Countries with four data points correspond to surveys collected in 2008, 2010, 2012, and 2014 (Argentina, Bolivia, Brazil and Paraguay), countries with five data points to 2006, 2008, 2010, 2012, and 2014 (Chile, Dominican Republic, Haiti, Peru, Uruguay and Venezuela), and countries with six data points correspond to surveys collected in 2004, 2006, 2008, 2010, 2012, and 2014 (Colombia, Costa Rica, Ecuador, El Salvador, Guatemala, Honduras, Mexico, Nicaragua, and Panama). 

plot_term_lrc_dt <- 
plot_term_lrc_dt %>%
  dplyr::select(pais, term, term_years_since_lrc, term_years_since_lrc_asinh, gdp_log, polity2, partipdem, delibdem, polyarchy)%>% 
  group_by(pais)%>%
  dplyr::summarize(term = mean(term, na.rm = T), term_years_since_lrc= mean(term_years_since_lrc,na.rm = T), term_years_since_lrc_asinh= mean(term_years_since_lrc_asinh,na.rm = T), gdp_log = mean(gdp_log, na..rm = T), polity2 = mean(polity2,na.rm = T),partipdem= mean(partipdem, na.rm = T), delibdem= mean(delibdem, na.rm = T), polyarchy= mean(polyarchy, na.rm = T))%>% 
  ungroup() 

plot_term_lrc_dt %>%
  mutate(
    abv = 
      paste(toupper(substr(str_remove(pais, " "), start = 1, stop = 3)))
  ) %>%
  ggplot() +
  aes(
    x = term_years_since_lrc_asinh,
    y = term
  ) +
  geom_point(shape = 16) +
  geom_smooth(method = "lm") +
  geom_text_repel(aes(label=abv), size=2.5,
  box.padding = ) +
  scale_x_continuous(breaks = asinh(custom_breaks), labels = custom_breaks) +
  scale_y_continuous(limits = c(0.5, 1.1), breaks = seq(0,1.1,by=0.2)) +
  labs(
    x = "Years since rule changed",
    y = "Proportion of correct answers",
    title = "How long is the presidential term?",
    caption = "Source: AmericasBarometer (2004-2014)"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot2_LA

print(plot2_LA)

###########################
# No. of representatives #
##########################

# Note on number of data points by country - For the question about the number of representatives in the lower house, countries have two datapoints (2012,2014), except for Suriname (only 2014).

plot_reps_right <- 
plot_reps_right %>%
  dplyr::select(pais, reps_right, reps_years_since, reps_years_since_asinh, gdp_log, polity2,partipdem, delibdem, polyarchy)%>% 
  group_by(pais)%>%
  dplyr::summarize(reps_right = mean(reps_right, na.rm = T), reps_years_since = mean(reps_years_since, na.rm = T), reps_years_since_asinh = mean(reps_years_since_asinh, na.rm = T), gdp_log = mean(gdp_log,na.rm = T), polity2 = mean(polity2,na.rm = T), partipdem= mean(partipdem, na.rm = T), delibdem= mean(delibdem, na.rm = T), polyarchy= mean(polyarchy, na.rm = T)) %>% 
  ungroup()

# Define the custom breaks
custom_breaks2 <- c(0, 15, 30, 45, 60, 75, 90)

plot_reps_right %>%
  mutate(
    abv = 
      paste(toupper(substr(str_remove(pais, " "), start = 1, stop = 3)))
  ) %>%
  ggplot() +
  aes(
    x = reps_years_since_asinh,
    y = reps_right
  ) +
  geom_point(shape = 16) +
  geom_smooth(method = "lm") +
  geom_text_repel(aes(label=abv), size=2.5,
  box.padding = ) +
  scale_x_continuous(breaks = asinh(custom_breaks2), labels = custom_breaks2) +
  scale_y_continuous(limits = c(-0.2, 0.7), breaks = seq(0,0.6,by=0.2)) +
  labs(
    x = "Years since number of representatives changed",
    y = "Proportion of correct answers",
    title = "How many representatives does the Lower House have?",
    caption = "Source: AmericasBarometer (2004-2014)"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot3_LA

print(plot3_LA)

####################################
##### Reelection term w/Lesotho ####
####################################

reelection_dt <- 
reelection_dt %>%
  dplyr::select(country, term, term_years_asinh, gdp_log, polity2, partipdem, delibdem, polyarchy)%>% 
  group_by(country)%>%
  dplyr::summarize(term = mean(term, na.rm = T), term_years_asinh = mean(term_years_asinh, na.rm = T), gdp_log = mean(gdp_log,na.rm = T), polity2 = mean(polity2,na.rm = T), partipdem= mean(partipdem, na.rm = T), delibdem= mean(delibdem, na.rm = T), polyarchy= mean(polyarchy, na.rm = T)) %>%
  ungroup() 

# Define the custom breaks 
custom_breaks3 <- c(0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25) 

reelection_dt%>%
  mutate(
    abv = 
      paste(toupper(substr(str_remove(country, " "), start = 1, stop = 3)))
  ) %>%
  ggplot() +
  aes(
    x = term_years_asinh,
    y = term
  ) +
  geom_point(shape = 16) +
  geom_smooth(method = "lm") +
  geom_text_repel(aes(label=abv), size=2.5,
  box.padding = ) + 
  scale_x_continuous(breaks = asinh(custom_breaks3), labels = custom_breaks3) +
  scale_y_continuous(breaks = seq(0,1,by=.2)) +
  labs(
    x = "Years since reelection rule changed",
    y = "Proportion of correct answers",
    title = "How many times someone can be elected President/PM?",
    caption= "Source: AfroBarometer (1999-2001)"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot1a_A

print(plot1a_A)


########################################
##### Reelection term w/out Lesotho ####
########################################

reelection_dt_nooutlier <- 
 reelection_dt %>%
 filter(country != "Lesotho") %>% # Apply the filter here
  dplyr::select(country, term, term_years_asinh, gdp_log, polity2, partipdem, delibdem, polyarchy)%>% 
  group_by(country)%>%
  dplyr::summarize(term = mean(term, na.rm = T), term_years_asinh = mean(term_years_asinh, na.rm = T), gdp_log = mean(gdp_log,na.rm = T), polity2 = mean(polity2,na.rm = T),  partipdem= mean(partipdem, na.rm = T), delibdem= mean(delibdem, na.rm = T), polyarchy= mean(polyarchy, na.rm = T)) %>% 
  ungroup() 

# Define the custom breaks 
custom_breaks4 <- c(0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5) 

reelection_dt_nooutlier %>%
  mutate(
    abv = 
      paste(toupper(substr(str_remove(country, " "), start = 1, stop = 3)))
  ) %>%
  ggplot() +
  aes(
    x = term_years_asinh,
    y = term
  ) +
  geom_point(shape = 16) +
  geom_smooth(method = "lm") +
  geom_text_repel(aes(label=abv), size=2.5,
  box.padding = ) + 
  scale_x_continuous(breaks = asinh(custom_breaks4), labels = custom_breaks4) +
  scale_y_continuous(limits = c(0.1, 0.8), breaks = seq(0.2,0.8,by=0.2)) +
  labs(
    x = "Years since reelection rule changed",
    y = "Proportion of correct answers",
    title = "How many times someone can be elected President/PM?",
    caption = "Source: Afrobarometer (1999-2001). Excludes Lesotho (outlier)"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot1b_A

print(plot1b_A)

# Arrange plots

common_theme <- theme(
  text = element_text(size = 13),  
  plot.caption = element_text(hjust = 1, size = 10), 
  plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
  axis.title = element_text(size = 14),  
  axis.text = element_text(size = 12),  
  legend.title = element_text(size = 10), 
  legend.text = element_text(size = 12)  
)

# Arrange plots
plot_political_institutions <- (plot1_LA | plot2_LA | plot3_LA | plot1b_A) + 
  plot_annotation(caption = "") & 
  common_theme

combined_plot <- plot_political_institutions + 
  plot_layout(ncol = 2, nrow = 2, guides = 'collect') &  
  ylab("Proportion of correct answers")

# Save the plot
ggsave(
  filename = "Figure3.tiff",
  plot = combined_plot,
  units = "in",
  height = 10,  
  width = 15  
)
```

```{r}
######################
#### Figure 4 ########
######################

#####################
###  Plur. party ####
#####################
party_dt <- 
party_dt %>%
  dplyr::select(country,party, party_years_asinh, gdp_log, polity2, partipdem, delibdem, polyarchy)%>% 
  group_by(country)%>%
  dplyr::summarize(party = mean(party, na.rm = T), party_years_asinh = mean(party_years_asinh, na.rm = T), gdp_log = mean(gdp_log,na.rm = T), polity2 = mean(polity2,na.rm = T), partipdem= mean(partipdem, na.rm = T), delibdem= mean(delibdem, na.rm = T), polyarchy= mean(polyarchy, na.rm = T)) %>%
  ungroup() 

# Define the custom breaks 
custom_breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40)

party_dt %>%
    mutate(
    abv = 
      paste(toupper(substr(str_remove(country, " "), start = 1, stop = 3)))
  ) %>%
  ggplot() +
  aes(
    x = party_years_asinh,
    y = party) +
  geom_point(shape = 16) +
  geom_smooth(method = "lm") +
  geom_text_repel(aes(label=abv), size=2.5,
  box.padding = ) + 
  scale_x_continuous(breaks = asinh(custom_breaks), labels = custom_breaks) +
  scale_y_continuous(breaks = seq(0,1,by=.2)) +
  labs(
    x = "Years since party control changed",
    y = "Proportion of correct answers",
    title = "Which party has the most seats in parliament/national assembly?",
    caption= "Source: Afrobarometer (1999-2001)"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot2_A

print(plot2_A)

#####################
## Vice president ###
#####################

leader_dt <- 
leader_dt %>%
  dplyr::select(country, leader, leader_years, gdp_log, polity2, partipdem, delibdem, polyarchy)%>% 
  group_by(country)%>%
  dplyr::summarize(leader = mean(leader, na.rm = T), leader_years= mean(leader_years,na.rm = T), gdp_log = mean(gdp_log,na.rm = T), polity2 = mean(polity2,na.rm = T), partipdem= mean(partipdem, na.rm = T), delibdem= mean(delibdem, na.rm = T), polyarchy= mean(polyarchy, na.rm = T)) %>%
  ungroup() 

leader_dt %>%
    mutate(
    abv = 
      paste(toupper(substr(str_remove(country, " "), start = 1, stop = 3)))
  ) %>%
  ggplot() +
  aes(
    x = leader_years,
    y = leader,
  ) +
  geom_point(shape = 16) +
  geom_smooth(method = "lm") +
  geom_text_repel(aes(label=abv), size=2.5,
  box.padding = ) + 
  scale_x_continuous(breaks = seq(0,15,by=2)) +
  scale_y_continuous(breaks = seq(0,1,by=.2)) +
  labs(
    x = "Years since deputy president/VP changed",
    y = "Proportion of correct answers",
    title = "What is the name of the deputy president/VP?",
    caption = "Source: Afrobarometer (1999-2001; 2005-2006)"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot3_A

print(plot3_A)

# Arrange plots
common_theme <- theme(
  text = element_text(size = 13),  
  plot.caption = element_text(hjust = 1, size = 10),  
  plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
  axis.title = element_text(size = 14), 
  axis.text = element_text(size = 12), 
  legend.title = element_text(size = 10),  
  legend.text = element_text(size = 12)  
)

# Arrange plots
plot_political_actors <- (plot2_A + plot3_A) + 
  plot_annotation(caption = "") & 
  common_theme

combined_plot_actors <- plot_political_actors + 
  plot_layout(ncol = 1, nrow = 2, guides = 'collect') &  
  ylab("Proportion of correct answers")

# Save the plot 
ggsave(
  filename = "Figure4.tiff",
  plot = combined_plot_actors,
  units = "in",
  height = 10,  
  width = 10  
)
```

```{r}
######################
#### Figure 5 ########
######################

#########################
###  Prog. structure ####
#########################

# Plot
cses_party_mean[which(cses_party_mean$DK_percentage <0.9),]%>%
  ggplot()+
  aes(
    x = prog_struc,
    y = dk_dummy
  ) +
  geom_point(shape=16) +
  geom_smooth(method = "lm") +
  scale_x_continuous(breaks = seq(0,1,by=.1)) +
  scale_y_continuous(breaks = seq(0,1,by=.2)) +
  labs(
    x = "Party programmatic structure",
    y = "Proportion of DK answers",
    title = "Where would you place party [X] on a L-R scale?"
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot1_CSES

print (plot1_CSES)

#########################
######## Party age ######
#########################

# Define the custom breaks 
custom_breaks <- c(0, 15, 30, 60, 75, 90, 120, 150, 180) 

cses_party_mean[which(cses_party_mean$DK_percentage <0.9),]%>%
  ggplot()+
  aes(
    x = party_age_asinh,
    y = dk_dummy
  ) +
  geom_point(shape=16) +
  geom_smooth(method = "lm") +
  scale_x_continuous(breaks = asinh(custom_breaks), labels = custom_breaks) +
  scale_y_continuous(breaks = seq(0,1,by=.2)) +
  labs(
    x = "Party age",
    y = "Proportion of DK answers",
  ) +
 theme_minimal()+
 theme(
    plot.title = element_text(face = "bold", hjust = 0, size= 14),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold", size=11, hjust = 1),
    axis.text.y = element_text(face = "bold", size=11),
    plot.caption = element_text(face = "italic", hjust = 1, size=9)
  ) -> plot2_CSES

print (plot2_CSES)


common_theme <- theme(
  text = element_text(size = 13),  
  plot.caption = element_text(face = "italic", hjust = 1, size = 10), 
  plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm"),
  axis.title = element_text(size = 14), 
  axis.text = element_text(size = 12),  
  legend.title = element_text(size = 10),  
  legend.text = element_text(size = 12)
  )

# Arrange plots 
combined_CSES_plot <- (plot1_CSES + plot2_CSES) + 
  plot_annotation(caption = "For party age data with x-axis relabeled using asinh function. Source: CSES data") & 
  common_theme

combined_CSES_plot  <- combined_CSES_plot  + 
  plot_layout(ncol = 1, nrow = 2, guides = 'collect') &  
  ylab("Proportion of DK answers")

# Save the plot 
ggsave(
  filename = "Figure5.tiff",
  plot = combined_CSES_plot,
  units = "in",
  height = 10,  
  width = 10   
)
```


```{r}
################################
#### Supplementary material ####
################################

################################
### Bivariate OLS regression ###
################################

##############################
### S4 - Number of states  ###
##############################

states_reg = lm(states ~ states_years_since_asinh, data = plot_states_dt)
# Robust standard errors
coeftest(states_reg, vcov. = vcovHC(states_reg, type = "HC0"))

# Output in HTML
stargazer(states_reg, type = "html", out="./SI/SI_Table4.html",
          se = list(sqrt(diag(vcovHC(states_reg, type = "HC0")))),
          title = "Bivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "States: Proportion Correct", 
          covariate.labels = c("Years Since (asinh)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

# Interpretation using asinh transformation
# Predicted values at the mean of the asinh-transformed variable and at one standard deviation above the mean (asinh_X_mean + asinh_SD):

mean(plot_states_dt$states_years_since_asinh) 
mean(plot_states_dt$states_years_since_asinh) + sd(plot_states_dt$states_years_since_asinh)
  
# Data from the bivariate regression model using robust SE
# Predicted_Y_mean = β₀ + β₁ * asinh_X_mean
0.344779 + 0.083256*3.785043 
# Predicted_Y_1SDabove = β₀ + β₁ * (asinh_X_mean + asinh_SD)
0.344779 + 0.083256*4.913791
# Difference between predicted values of Y
0.7538816-0.6599065

##############################
### S5 - Presidential term ###
##############################

term_lrc_reg = lm(term ~ term_years_since_lrc_asinh, data = plot_term_lrc_dt)
# Robust standard errors
coeftest(term_lrc_reg, vcov. = vcovHC(term_lrc_reg, type = "HC0"))
# Output in HTML
stargazer(term_lrc_reg, type = "html", out="./SI/SI_Table5.html",
          se = list(sqrt(diag(vcovHC(term_lrc_reg, type = "HC0")))),
          title = "Bivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Term Length: Proportion Correct", 
          covariate.labels = c("Years Since (asinh)","Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

# Interpretation using asinh transformation
# Predicted values at the mean of the asinh-transformed variable and at one standard deviation above the mean (asinh_X_mean + asinh_SD):

mean(plot_term_lrc_dt$term_years_since_lrc_asinh)
mean(plot_term_lrc_dt$term_years_since_lrc_asinh) + sd(plot_term_lrc_dt$term_years_since_lrc_asinh)

# Data from the bivariate regression model using robust SE
# Predicted_Y_mean = β₀ + β₁ * asinh_X_mean
0.70845 + 0.03792*4.178305
# Predicted_Y_1SDabove = β₀ + β₁ * (asinh_X_mean + asinh_SD)
0.70845 + 0.03792*5.136863
# Difference between predicted values of Y
0.9032398-0.8668913

################################
#  S6 - No. of representatives #
################################

reps_reg = lm(reps_right ~ reps_years_since_asinh, data = plot_reps_right)
# Robust standard errors
coeftest(reps_reg, vcov. = vcovHC(reps_reg, type = "HC0"))
# Output in HTML
stargazer(reps_reg, type = "html", out="./SI/SI_Table6.html",
          se = list(sqrt(diag(vcovHC(reps_reg, type = "HC0")))),
          title = "Bivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Representatives: Proportion Correct",           
          covariate.labels = c("Years Since (asinh)","Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

# Interpretation using asinh transformation
# Predicted values at the mean of the asinh-transformed variable and at one standard deviation above the mean (asinh_X_mean + asinh_SD):

mean(plot_reps_right$reps_years_since_asinh) 
mean(plot_reps_right$reps_years_since_asinh) + sd(plot_reps_right$reps_years_since_asinh)

# Data from the bivariate regression model using robust SE
# Predicted_Y_mean = β₀ + β₁ * asinh_X_mean
0.042882+ 0.058546*3.012066
# Predicted_Y_1SDabove = β₀ + β₁ * (asinh_X_mean + asinh_SD)
0.042882 + 0.058546*4.366207
# Difference between predicted values of Y
0.298506 - 0.2192264

##############################################
##### S7 - Reelection term w/out Lesotho ####
#############################################

reelection_reg_alt <- lm(term ~ term_years_asinh, data = reelection_dt, subset = (country != "Lesotho"))
# Robust standard errors
coeftest(reelection_reg_alt, vcov. = vcovHC(reelection_reg_alt, type = "HC0"))
# Output in HTML
stargazer(reelection_reg_alt, type = "html", out="./SI/SI_Table7.html",
          se = list(sqrt(diag(vcovHC(reelection_reg_alt, type = "HC0")))),
          title = "Bivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Terms: Proportion Correct", 
          covariate.labels = c("Years Since (asinh)","Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes")))

# Interpretation using asinh transformation
# Predicted values at the mean of the asinh-transformed variable and at one standard deviation above the mean (asinh_X_mean + asinh_SD):

mean(reelection_dt$term_years_asinh[reelection_dt$country != "Lesotho"], na.rm = TRUE)
mean(reelection_dt$term_years_asinh[reelection_dt$country != "Lesotho"], na.rm = TRUE)+ sd(reelection_dt$term_years_asinh[reelection_dt$country != "Lesotho"], na.rm = TRUE)
  
# Data from the bivariate regression model using robust SE
# Predicted_Y_mean = β₀ + β₁ * asinh_X_mean
0.262732 + 0.057471*3.135457
# Predicted_Y_1SDabove = β₀ + β₁ * (asinh_X_mean + asinh_SD)
0.262732  + 0.057471*3.537231
# Difference between predicted values of Y
0.4660202 - 0.4429298

##########################
###  S8 - Plur. party ####
##########################

plurality_reg = lm(party ~ party_years_asinh, data = party_dt)
#Robust standard errors
coeftest(plurality_reg, vcov. = vcovHC(plurality_reg, type = "HC0"))
# Output in HTML
stargazer(plurality_reg, type = "html", out="./SI/SI_Table8.html",
          se = list(sqrt(diag(vcovHC(plurality_reg, type = "HC0")))),
          title = "Bivariate regression",
          header = FALSE, 
          model.names = FALSE,
          dep.var.labels = "Party: Proportion Correct",
          covariate.labels = c("Years Since (asinh)","Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

# Interpretation using asinh transformation
# Predicted values at the mean of the asinh-transformed variable and at one standard deviation above the mean (asinh_X_mean + asinh_SD):

mean(party_dt$party_years_asinh)
mean(party_dt$party_years_asinh) + sd(party_dt$party_years_asinh)
  
# Data from the bivariate regression model using robust SE
# Predicted_Y_mean = β₀ + β₁ * asinh_X_mean
0.211042 + 0.158096*2.643404
# Predicted_Y_1SDabove = β₀ + β₁ * (asinh_X_mean + asinh_SD)
0.211042 + 0.158096*3.693714
# Difference between predicted values of Y
0.7950034-0.6289536

#############################
###  S9 - Vice president ####
#############################
leader_reg = lm(leader ~ leader_years, data = leader_dt)

#Robust standard errors
coeftest(leader_reg, vcov. = vcovHC(leader_reg, type = "HC0"))

# Output in HTML
stargazer(leader_reg, type = "html", out="./SI/SI_Table9.html",
          se = list(sqrt(diag(vcovHC(leader_reg, type = "HC0")))),
          title = "Bivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Vice President: Proportion Correct",
          covariate.labels = c("Years Since","Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

########################################
### S10 - Reelection term w/ Lesotho ###
########################################

reelection_reg = lm(term ~ term_years_asinh, data = reelection_dt)
summary(reelection_reg)
# Robust standard errors
coeftest(reelection_reg, vcov. = vcovHC(reelection_reg, type = "HC0"))
# Output in HTML
stargazer(reelection_reg, type = "html", out="./SI/SI_Table10.html",
          se = list(sqrt(diag(vcovHC(reelection_reg, type = "HC0")))),
          title = "Bivariate regression",
          header = FALSE,
          model.names = FALSE, 
          dep.var.labels = "Terms: Proportion Correct",
          covariate.labels = c("Years Since (asinh)","Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 
```


```{r}
######################################
  ### Multivariate OLS regression ###
#####################################

##############################
### S11 - Number of states ###
##############################

states_reg_alt = lm(states ~ states_years_since_asinh + gdp_log + polity2, data = plot_states_dt)
# Robust standard errors
coeftest(states_reg_alt, vcov. = vcovHC(states_reg_alt, type = "HC0"))
# Output in HTML
stargazer(states_reg_alt, type = "html", out="./SI/SI_Table11.html",
          se = list(sqrt(diag(vcovHC(states_reg_alt, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "States: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

############################################
######## VDEM controls (S11A - S11C) #######
############################################

states_reg_1 = lm(states ~ states_years_since_asinh + gdp_log + partipdem, data = plot_states_dt)
# robust stardard errors
coeftest(states_reg_1, vcov. = vcovHC(states_reg_1, type = "HC0"))
# Output in HTML
stargazer(states_reg_1, type = "html", out="./SI/SI_Table11A.html",
          se = list(sqrt(diag(vcovHC(states_reg_1, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "States: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Participation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

states_reg_2 = lm(states ~ states_years_since_asinh + gdp_log + delibdem, data = plot_states_dt)
#robust standard errors
coeftest(states_reg_2, vcov. = vcovHC(states_reg_2, type = "HC0"))
# Output in HTML
stargazer(states_reg_2, type = "html", out="./SI/SI_Table11B.html",
          se = list(sqrt(diag(vcovHC(states_reg_2, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "States: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Deliberation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes")))

states_reg_3 = lm(states ~ states_years_since_asinh + gdp_log + polyarchy, data = plot_states_dt)
# Robust standard errors
coeftest(states_reg_3, vcov. = vcovHC(states_reg_3, type = "HC0"))
# Output in HTML
stargazer(states_reg_3, type = "html", out="./SI/SI_Table11C.html",
          se = list(sqrt(diag(vcovHC(states_reg_3, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "States: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polyarchy(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

###############################
### S12 - Presidential term ###
###############################
term_lrc_reg_alt = lm(term ~ term_years_since_lrc_asinh + gdp_log + polity2, data = plot_term_lrc_dt)
#Robust standard errors
coeftest(term_lrc_reg_alt, vcov. = vcovHC(term_lrc_reg_alt, type = "HC0"))
# Output in HTML
stargazer(term_lrc_reg_alt, type = "html", out="./SI/SI_Table12.html",
          se = list(sqrt(diag(vcovHC(term_lrc_reg_alt, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Term Length: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

############################################
######## VDEM controls (S12A - S12C) #######
############################################

term_reg_1 = lm(term ~ term_years_since_lrc_asinh + gdp_log + partipdem, data = plot_term_lrc_dt)
#Robust standard errors
coeftest(term_reg_1, vcov. = vcovHC(term_reg_1, type = "HC0"))
# Output in HTML
stargazer(term_reg_1, type = "html", out="./SI/SI_Table12A.html",
          se = list(sqrt(diag(vcovHC(term_reg_1, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Term Length: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Participation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes")))


term_reg_2 = lm(term ~term_years_since_lrc_asinh + gdp_log + delibdem, data = plot_term_lrc_dt)
#Robust standard error
coeftest(term_reg_2, vcov. = vcovHC(term_reg_2, type = "HC0"))
# Output in HTML
stargazer(term_reg_2, type = "html", out="./SI/SI_Table12B.html",
          se = list(sqrt(diag(vcovHC(term_reg_2, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Term Length: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Deliberation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 


term_reg_3 = lm(term ~ term_years_since_lrc_asinh + gdp_log + polyarchy, data = plot_term_lrc_dt)
#Robust standard errors
coeftest(term_reg_3, vcov. = vcovHC(term_reg_3, type = "HC0"))
# Output in HTML
stargazer(term_reg_3, type = "html", out="./SI/SI_Table12C.html",
          se = list(sqrt(diag(vcovHC(term_reg_3, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE,
          dep.var.labels = "Term Length: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polyarchy(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

####################################
### S13 - No. of representatives ###
####################################

reps_reg_alt= lm(reps_right ~ reps_years_since_asinh + gdp_log + polity2, data = plot_reps_right)
#Robust standard errors
coeftest(reps_reg_alt, vcov. = vcovHC(reps_reg_alt, type = "HC0"))
# Output in HTML
stargazer(reps_reg_alt, type = "html", out="./SI/SI_Table13.html",
          se = list(sqrt(diag(vcovHC(reps_reg_alt, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Representatives: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

############################################
######## VDEM controls (S13A - S13C) #######
############################################

reps_reg_1 = lm(reps_right ~ reps_years_since_asinh + gdp_log + partipdem, data = plot_reps_right)
# Robust standard errors
coeftest(reps_reg_1, vcov. = vcovHC(reps_reg_1, type = "HC0"))
# Output in HTML
stargazer(reps_reg_1, type = "html", out="./SI/SI_Table13A.html",
          se = list(sqrt(diag(vcovHC(reps_reg_1, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Representatives: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Participation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

reps_reg_2 = lm(reps_right ~ reps_years_since_asinh + gdp_log + delibdem, data = plot_reps_right)
# Robust standard errors
coeftest(reps_reg_2, vcov. = vcovHC(reps_reg_2, type = "HC0"))
# Output in HTML
stargazer(reps_reg_2, type = "html", out="./SI/SI_Table13B.html",
          se = list(sqrt(diag(vcovHC(reps_reg_2, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Representatives: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Deliberation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

reps_reg_3 = lm(reps_right ~ reps_years_since_asinh + gdp_log + polyarchy, data = plot_reps_right)
# Robust standard errors
coeftest(reps_reg_3, vcov. = vcovHC(reps_reg_3, type = "HC0"))
# Output in HTML
stargazer(reps_reg_3, type = "html", out="./SI/SI_Table13C.html",
          se = list(sqrt(diag(vcovHC(reps_reg_3, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE,
          dep.var.labels = "Representatives: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polyarchy(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

###########################################
### S14 - Reelection term w/out Lesotho ###
###########################################

reelection_reg_alt = lm(term ~ term_years_asinh + gdp_log + polity2, data = reelection_dt, subset = (country != "Lesotho"))
#Robust standard errors
coeftest(reelection_reg_alt, vcov. = vcovHC(reelection_reg_alt, type = "HC0"))
# Output in HTML
stargazer(reelection_reg_alt,type = "html", out="./SI/SI_Table14.html",
          se = list(sqrt(diag(vcovHC(reelection_reg_alt, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE,
          dep.var.labels = "Terms: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

############################################
######## VDEM controls (S14A - S14C) #######
############################################

reelection_reg_alt1 = lm(term ~ term_years_asinh + gdp_log + partipdem, data = reelection_dt, subset = (country != "Lesotho"))
# Robust standard errors
coeftest(reelection_reg_alt1, vcov. = vcovHC(reelection_reg_alt1, type = "HC0"))
# Output in HTML
stargazer(reelection_reg_alt1,type = "html", out="./SI/SI_Table14A.html",
          se = list(sqrt(diag(vcovHC(reelection_reg_alt1, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Terms: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Participation (VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 


reelection_reg_alt2 = lm(term ~ term_years_asinh + gdp_log + delibdem, data = reelection_dt, subset = (country != "Lesotho"))
# Robust standard errors
coeftest(reelection_reg_alt2, vcov. = vcovHC(reelection_reg_alt2, type = "HC0"))
# Output in HTML
stargazer(reelection_reg_alt2,type = "html", out="./SI/SI_Table14B.html",
          se = list(sqrt(diag(vcovHC(reelection_reg_alt2, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Terms: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Deliberation (VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

reelection_reg_alt3 = lm(term ~ term_years_asinh + gdp_log + polyarchy, data = reelection_dt, subset = (country != "Lesotho"))
#Robust standard errors
coeftest(reelection_reg_alt3, vcov. = vcovHC(reelection_reg_alt3, type = "HC0"))
# Output in HTML
stargazer(reelection_reg_alt3,type = "html", out="./SI/SI_Table14C.html",
          se = list(sqrt(diag(vcovHC(reelection_reg_alt3, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE,
          dep.var.labels = "Terms: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polyarchy (VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

##########################
###  S15 - Plur.party ####
##########################

plurality_reg_alt = lm(party ~ party_years_asinh + gdp_log + polity2, data = party_dt)
#Robust standard errors
coeftest(plurality_reg_alt, vcov. = vcovHC(plurality_reg_alt, type = "HC0"))
# Output in HTML
stargazer(plurality_reg_alt, type = "html", out="./SI/SI_Table15.html",
          se = list(sqrt(diag(vcovHC(plurality_reg_alt, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Party: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 


############################################
######## VDEM controls (S15A - S15C) #######
############################################

plurality_reg_1 = lm(party ~ party_years_asinh + gdp_log + partipdem, data = party_dt)
# Robust standard errors
coeftest(plurality_reg_1, vcov. = vcovHC(plurality_reg_1, type = "HC0"))
# Output in HTML
stargazer(plurality_reg_1, type = "html", out="./SI/SI_Table15A.html",
          se = list(sqrt(diag(vcovHC(plurality_reg_1, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Party: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Participation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

plurality_reg_2 = lm(party ~ party_years_asinh + gdp_log + delibdem, data = party_dt)
#Robust standard errors
coeftest(plurality_reg_2, vcov. = vcovHC(plurality_reg_2, type = "HC0"))
# Output in HTML
stargazer(plurality_reg_2, type = "html", out="./SI/SI_Table15B.html",
          se = list(sqrt(diag(vcovHC(plurality_reg_2, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Party: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Deliberation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 


plurality_reg_3 = lm(party ~ party_years_asinh + gdp_log + polyarchy, data = party_dt)
# Robust standard errors
coeftest(plurality_reg_3, vcov. = vcovHC(plurality_reg_3, type = "HC0"))
# Output in HTML
stargazer(plurality_reg_3, type = "html", out="./SI/SI_Table15C.html",
          se = list(sqrt(diag(vcovHC(plurality_reg_3, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Party: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polyarchy (VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

##############################
###  S16 - Vice president ####
##############################

leader_reg_alt = lm(leader ~ leader_years + gdp_log + polity2, data = leader_dt)
#Robust standard errors
coeftest(leader_reg_alt, vcov. = vcovHC(leader_reg_alt, type = "HC0"))
# Output in HTML
stargazer(leader_reg_alt, type = "html", out="./SI/SI_Table16.html",
          se = list(sqrt(diag(vcovHC(leader_reg_alt, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE,
          model.names = FALSE, 
          dep.var.labels = "Vice President: Proportion Correct",
          covariate.labels = c("Years Since", "GDP per capita (log)", "Polity V", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

############################################
######## VDEM controls (S16A - S16C) #######
############################################

leader_reg_1 = lm(leader ~ leader_years + gdp_log + partipdem, data = leader_dt)
# Robust standard errors
coeftest(leader_reg_1, vcov. = vcovHC(leader_reg_1, type = "HC0"))
# Output in HTML
stargazer(leader_reg_1, type = "html", out="./SI/SI_Table16A.html",
          se = list(sqrt(diag(vcovHC(leader_reg_1, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Vice President: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Participation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

leader_reg_2 = lm(leader ~ leader_years + gdp_log + delibdem, data = leader_dt)
#Robust standard errors
coeftest(leader_reg_2, vcov. = vcovHC(leader_reg_2, type = "HC0"))
# Output in HTML
stargazer(leader_reg_2, type = "html", out="./SI/SI_Table16B.html",
          se = list(sqrt(diag(vcovHC(leader_reg_2, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Vice President: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Deliberation(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 

leader_reg_3 = lm(leader ~ leader_years + gdp_log + polyarchy, data = leader_dt)
#Robust standard errors
coeftest(leader_reg_3, vcov. = vcovHC(leader_reg_3, type = "HC0"))
# Output in HTML
stargazer(leader_reg_3, type = "html", out="./SI/SI_Table16C.html",
          se = list(sqrt(diag(vcovHC(leader_reg_3, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Vice President: Proportion Correct",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polyarchy(VDEM)", "Intercept"),
          add.lines = list(c("Robust Standard Errors", "Yes"))) 
```

```{r}
####################################################
####### S17 - Party programmatic structure #########
####################################################

m3 = lm(dk_dummy ~ prog_struc + partysize + survey + mail_any + internet_any + phone_any + person_any, data = cses_party_mean[cses_party_mean$DK_percentage <0.9, ])
summary(m3)
# Robust standard errors
coeftest(m3, vcov. = vcovHC(m3, type = "HC0"))
# Output in HTML
stargazer(m3, type = "html", out="./SI/SI_Table17.html",
          se = list(sqrt(diag(vcovHC(m3, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Proportion “Don’t Know” Answers",
          covariate.labels = c("Prog structure", "Party size"),
          add.lines = list(c("Robust Standard Errors", "Yes")),
          scale_down = 0.3  
) 


#################################
####### S18 - Party age #########
#################################

m5 = lm(dk_dummy ~ party_age_asinh + partysize + survey + mail_any + internet_any + phone_any + person_any, data = cses_party_mean[cses_party_mean$DK_percentage <0.9, ])
# Robust standard errors
coeftest(m5, vcov. = vcovHC(m5, type = "HC0"))
# Output in HTML
stargazer(m5, type = "html", out="./SI/SI_Table18.html",
          se = list(sqrt(diag(vcovHC(m5, type = "HC0")))),
          title = "Multivariate regression",
          header = FALSE, 
          model.names = FALSE, 
          dep.var.labels = "Proportion “Don’t Know” Answers",
          covariate.labels = c("Party age (asinh)", "Party size"),
          add.lines = list(c("Robust Standard Errors", "Yes")),
          scale_down = 0.3  
) 

# Coefficients with asinh transformations 
mean.age <- mean(cses_party_mean$party_age_asinh, na.rm=T)
coef(m5)["party_age_asinh"]/sqrt(mean.age^2 + 1)
sd(cses_party_mean$party_age, na.rm=T)*coef(m5)["party_age_asinh"]/sqrt(mean.age^2 + 1)
```

```{r}
#########################################
###### Hierarchical Linear Models #######
#########################################

###########################################
###### S19A - S19B.Number of states  ######
###########################################

# Random Intercept
hlm_states <- lmer(states ~ states_years_since_asinh + gdp_log + polity2 + age + gender + edu +
                       (1 | pais/countrywave),
                       data = dt2)

stargazer(hlm_states, type = "html", out="./SI/SI_Table19A.html",
          title = "Hierarchical Linear Model Results using Random Intercept",
          dep.var.labels = "States: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

#Random slopes
hlm_states_alt <- lmer(
  states ~ states_years_since_asinh + gdp_log + polity2 + age + gender + edu + 
    (1 + age + gender + edu | pais/countrywave), 
  data = dt2
)

#Output in HTML
stargazer(hlm_states_alt, type = "html", out="./SI/SI_Table19B.html",
          title = "Hierarchical Linear Model Results using Random Slopes",
          dep.var.labels = "States: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

####################################
### S20A-S20B. Presidential term ###
####################################

# Random Intercept
hlm_term<- lmer(term ~ term_years_since_lrc_asinh + gdp_log + polity2 + age + gender + edu +
                       (1 | pais/countrywave),
                       data = dt2)
# Output in HTML
stargazer(hlm_term, type = "html", out="./SI/SI_Table20A.html",
          title = "Hierarchical Linear Model Results using Random Intercept",
          dep.var.labels = "Term Length: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

# Random Slopes
hlm_term_alt <- lmer(
  term ~ age + gender + edu + term_years_since_lrc_asinh + gdp_log + polity2 +
    (1 + age + gender + edu | pais/countrywave), 
  data = dt2
)

# Output in HTML
stargazer(hlm_term_alt, type = "html", out="./SI/SI_Table20B.html",
          title = "Hierarchical Linear Model Results using Random Slopes",
          dep.var.labels = "Term Length: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

#####################################
#  S21A-21B. No. of representatives #
#####################################

# Random Intercept
hlm_reps <- lmer(reps_right ~ reps_years_since_asinh + gdp_log + polity2 + age + gender + edu + 
                       (1 | pais/countrywave),
                       data = dt2)
# Output in HTML
stargazer(hlm_reps, type = "html", out="./SI/SI_Table21A.html",
          title = "Hierarchical Linear Model Results using Random Intercept",
          dep.var.labels = "Representatives: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

# Random Slopes
hlm_reps_alt <- lmer(
 reps_right ~ reps_years_since_asinh + gdp_log + polity2 + age + gender + edu + 
    (1 + age + gender + edu | pais/countrywave), 
  data = dt2
)

# Output in HTML
stargazer(hlm_reps_alt, type = "html", out="./SI/SI_Table21B.html",
          title = "Hierarchical Linear Model Results using Random Slopes",
          dep.var.labels = "Representatives: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))


###################################################
##### S22A-22B - Reelection term w/out Lesotho ####
###################################################

#Random Intercept
hlm_reelection_nooutlier <- lmer(
  term ~ term_years_asinh + gdp_log + polity2 + age + gender + edu + (1 | country/countryyear),
  data = dt,
  subset = (country != "Lesotho")
)

# Output in HTML
stargazer(hlm_reelection_nooutlier, type = "html", out="./SI/SI_Table22A.html",
          title = "Hierarchical Linear Model Results using Random Intercept (excludes Lesotho)",
          dep.var.labels = "Terms: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

#Random Slopes
hlm_reelection_nooutlier_alt <- lmer(
  term ~ term_years_asinh + gdp_log + polity2 + age + gender + edu + 
    (1 + age + gender + edu | country/countryyear), 
  data = dt,
  subset = (country != "Lesotho")
)

# Output in HTML
stargazer(hlm_reelection_nooutlier, type = "html", out="./SI/SI_Table22B.html",
          title = "Hierarchical Linear Model Results using Random Slopes (excludes Lesotho)",
          dep.var.labels = "Terms: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

################################
###  S23A-23B - Plur. party ####
################################

# Random Intercept
hlm_plur_party <- lmer(party ~ party_years_asinh + gdp_log + polity2 + age + gender + edu + 
                       (1 | country/countryyear),
                       data = dt)
# Output in HTML
stargazer(hlm_plur_party, type = "html", out="./SI/SI_Table23A.html",
          title = "Hierarchical Linear Model Results using Random Intercept",
          dep.var.labels = "Plurality Party: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

# Random Slopes
hlm_plur_party_alt <- lmer(
  party ~ party_years_asinh + gdp_log + polity2 + age + gender + edu +
    (1 + age + gender + edu | country/countryyear), 
  data = dt
)

# Output in HTML
stargazer(hlm_plur_party_alt, type = "html", out="./SI/SI_Table23B.html",
          title = "Hierarchical Linear Model Results using Random Slopes",
          dep.var.labels = "Plurality Party: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))


###################################
###  S24A-24B - Vice president ####
###################################

# Random Intercept
hlm_leader <- lmer(leader ~ leader_years + gdp_log + polity2 + age + gender + edu +
                       (1 | country/countryyear),
                       data = dt)

# Output in HTML
stargazer(hlm_leader, type = "html", out="./SI/SI_Table24A.html",
          title = "Hierarchical Linear Model Results using Random Intercept",
          dep.var.labels = "Vice President: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

# Random Slopes
hlm_leader_alt <- lmer(
  leader ~ leader_years + gdp_log + polity2 + age + gender + edu +
    (1 + age + gender + edu | country/countryyear), 
  data = dt
)

# Output in HTML
stargazer(hlm_leader_alt, type = "html", out="./SI/SI_Table24B.html",
          title = "Hierarchical Linear Model Results using Random Slopes",
          dep.var.labels = "Vice President: Correct (0/1)",
          covariate.labels = c("Years Since (asinh)", "GDP per capita (log)", "Polity V", "Age", "Gender", "Education"))

####################################################
#### S25A-25B - Party programmatic structure #######
####################################################

# Random Intercept
hlm_prog_struc <- lmer(dk_dummy ~ prog_struc + partysize + age + gender + edu + (1 | survey/party),
              data = dt3[dt3$DK_percentage <0.9, ])

# Output in HTML
stargazer(hlm_prog_struc, type = "html", out="./SI/SI_Table25A.html",
          title = "Hierarchical Linear Model Results using Random Intercept",
          dep.var.labels = "Proportion “Don’t Know” Answers",
          covariate.labels = c("Party prog structure", "Party size", "Age", "Gender", "Education"))

# Random Slopes
hlm_prog_struc_alt <- lmer(
  dk_dummy ~ prog_struc + partysize + age + gender + edu + 
    (1 + age + gender + edu | survey/party), 
  data = dt3[dt3$DK_percentage <0.9, ]
)

# Output in HTML
stargazer(hlm_prog_struc_alt, type = "html", out="./SI/SI_Table25B.html",
          title = "Hierarchical Linear Model Results using Random Slopes",
          dep.var.labels = "Proportion “Don’t Know” Answers",
          covariate.labels = c("Party prog structure", "Party size", "Age", "Gender", "Education"))

#################################
##### S26A-26B - Party age ######
#################################

# Random Intercept
hlm_age <- lmer(dk_dummy ~ party_age_asinh + partysize + age + gender + edu + (1 | survey/party),
              data = dt3[dt3$DK_percentage <0.9, ])

# Output in HTML
stargazer(hlm_age, type = "html", out="./SI/SI_Table26A.html",
          title = "Hierarchical Linear Model Results using Random Intercept",
          dep.var.labels = "Proportion “Don’t Know” Answers",
          covariate.labels = c("Party Age (asinh)", "Party size", "Age", "Gender", "Education"))

# Random Slopes
hlm_age_alt <- lmer(
  dk_dummy ~ party_age_asinh + partysize + age + gender + edu + 
    (1 + age + gender + edu | survey/party), 
  data = dt3[dt3$DK_percentage <0.9, ]
)

# Output in HTML
stargazer(hlm_age_alt, type = "html", out="./SI/SI_Table26B.html",
          title = "Hierarchical Linear Model Results using Random Slopes",
          dep.var.labels = "Proportion “Don’t Know” Answers",
          covariate.labels = c("Party Age (asinh)", "Party size", "Age", "Gender", "Education"))
```

